Importing Data
wholesale_computer = scan("~/desktop/computerwholesale.csv")
wholesale_computer_timeseries = ts(wholesale_computer, frequency = 12, start = c(1997, 1))
Calling the packages needed
library(ggplot2)
library(tseries)
## Warning: package 'tseries' was built under R version 3.3.2
library(TSA)
## Loading required package: leaps
## Warning: package 'leaps' was built under R version 3.3.2
## Loading required package: locfit
## locfit 1.5-9.1 2013-03-22
## Loading required package: mgcv
## Loading required package: nlme
## This is mgcv 1.8-12. For overview type 'help("mgcv-package")'.
##
## Attaching package: 'TSA'
## The following objects are masked from 'package:stats':
##
## acf, arima
## The following object is masked from 'package:utils':
##
## tar
library(forecast)
## Warning: package 'forecast' was built under R version 3.3.2
## Warning in as.POSIXlt.POSIXct(Sys.time()): unknown timezone 'zone/tz/2018e.
## 1.0/zoneinfo/Asia/Shanghai'
##
## Attaching package: 'forecast'
## The following object is masked from 'package:nlme':
##
## getResponse
library(quantmod)
## Warning: package 'quantmod' was built under R version 3.3.2
## Loading required package: xts
## Warning: package 'xts' was built under R version 3.3.2
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 3.3.2
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
## Loading required package: TTR
## Warning: package 'TTR' was built under R version 3.3.2
## Version 0.4-0 included new data defaults. See ?getSymbols.
library(rugarch)
## Warning: package 'rugarch' was built under R version 3.3.2
## Loading required package: parallel
##
## Attaching package: 'rugarch'
## The following object is masked from 'package:stats':
##
## sigma
Initial Plot
plot.ts(wholesale_computer_timeseries)

Decompose data into trend, seasonal and remainder components
fit = stl(wholesale_computer_timeseries, s.window = "periodic")
plot(fit)

Spectral Analysis
raw_spec = spec.pgram(wholesale_computer_timeseries, taper = 0)
plot(raw_spec)

plot(raw_spec, log = "no")
abline(v = 1/7, col = "red")
abline(v = 7/24, col = "red")
abline(v = 9/24, col = "red")
text(x = c(0.2, 0.3, 0.4), labels = c("1/7", "7/24", "9/24"), y = c(1500000, 1700000, 2000000), col = "red")

sp = lm(formula = diff(wholesale_computer_timeseries) ~ cos(2*pi*1/7*seq(from = 1, to = 247)) + sin(2*pi*1/7*seq(from = 1, to = 247)) + cos(2*pi*7/24*seq(from = 1, to = 247)) + sin(2*pi*7/24*seq(from = 1, to = 247)) +cos(2*pi*9/24*seq(from = 1, to = 247)) + sin(2*pi*9/24*seq(from = 1, to = 247)), diff(wholesale_computer_timeseries))
curve(38.7 - 26.334*cos(2*pi*1/7*x) +10.967*sin(2*pi*1/7*x) + 11.275*cos(2*pi*7/24*x) -9.285*sin(2*pi*7/24*x) -2.8 *cos(2*pi*9/24*x) +68.481*sin(2*pi*9/24*x), from = 1, to = 247, ylab = "spectral")
lines(diff(wholesale_computer), col = "blue")

qqnorm(window(rstandard(sp)));qqline(window(rstandard(sp)))

shapiro.test(rstandard(sp))
##
## Shapiro-Wilk normality test
##
## data: rstandard(sp)
## W = 0.98288, p-value = 0.004535
ARIMA Model
#adf test
adf.test(wholesale_computer_timeseries, alternative="stationary", k=0)
##
## Augmented Dickey-Fuller Test
##
## data: wholesale_computer_timeseries
## Dickey-Fuller = -3.294, Lag order = 0, p-value = 0.07274
## alternative hypothesis: stationary
adf.test(diff(wholesale_computer_timeseries), alternative = "stationary", k=0)
## Warning in adf.test(diff(wholesale_computer_timeseries), alternative =
## "stationary", : p-value smaller than printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: diff(wholesale_computer_timeseries)
## Dickey-Fuller = -19.192, Lag order = 0, p-value = 0.01
## alternative hypothesis: stationary
acf(diff(wholesale_computer_timeseries))

pacf(diff(wholesale_computer_timeseries))

eacf(diff(wholesale_computer_timeseries))
## AR/MA
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 0 x x x o o x o o o x x o o o
## 1 x x x o o x o o o o x x o o
## 2 x x o o o o o o o o o x o o
## 3 o x x o o o o o o o o x o o
## 4 o o x o o o o o o o o x o o
## 5 x x x o o o o o o o o o o o
## 6 x o x x o o o o o o x o o o
## 7 x o o x o o o o o o o x o o
mean(diff(wholesale_computer_timeseries))
## [1] 38.74494
#plot to test for stationary
plot.ts(wholesale_computer_timeseries)

plot.ts(diff(wholesale_computer_timeseries))

#ARIMA model
arima(diff(wholesale_computer_timeseries), order = c(2, 1, 2), include.mean = FALSE)
##
## Call:
## arima(x = diff(wholesale_computer_timeseries), order = c(2, 1, 2), include.mean = FALSE)
##
## Coefficients:
## ar1 ar2 ma1 ma2
## -0.6591 -0.2978 -0.5741 -0.4259
## s.e. 0.1466 0.0628 0.1448 0.1444
##
## sigma^2 estimated as 193392: log likelihood = -1849.43, aic = 3706.86
arima(diff(wholesale_computer_timeseries), order = c(0, 1, 3), include.mean = FALSE)
##
## Call:
## arima(x = diff(wholesale_computer_timeseries), order = c(0, 1, 3), include.mean = FALSE)
##
## Coefficients:
## ma1 ma2 ma3
## -1.2243 0.1870 0.0373
## s.e. 0.0789 0.1462 0.0837
##
## sigma^2 estimated as 202130: log likelihood = -1854.79, aic = 3715.59
arima(diff(wholesale_computer_timeseries), order = c(1, 1, 3), include.mean = FALSE)
##
## Call:
## arima(x = diff(wholesale_computer_timeseries), order = c(1, 1, 3), include.mean = FALSE)
##
## Coefficients:
## ar1 ma1 ma2 ma3
## -0.5310 -0.6605 -0.5587 0.2192
## s.e. 0.1824 0.1756 0.1949 0.0564
##
## sigma^2 estimated as 198717: log likelihood = -1852.72, aic = 3713.44
arima(diff(wholesale_computer_timeseries), order = c(2, 1, 3), include.mean = FALSE)
##
## Call:
## arima(x = diff(wholesale_computer_timeseries), order = c(2, 1, 3), include.mean = FALSE)
##
## Coefficients:
## ar1 ar2 ma1 ma2 ma3
## -0.6641 -0.7899 -0.5453 0.1587 -0.6134
## s.e. 0.0999 0.0912 0.1240 0.1699 0.1274
##
## sigma^2 estimated as 186715: log likelihood = -1845.06, aic = 3700.12
arima(diff(wholesale_computer_timeseries), order = c(3, 1, 3), include.mean = FALSE)
##
## Call:
## arima(x = diff(wholesale_computer_timeseries), order = c(3, 1, 3), include.mean = FALSE)
##
## Coefficients:
## ar1 ar2 ar3 ma1 ma2 ma3
## -0.6574 -0.7851 0.0039 -0.5506 0.1621 -0.6115
## s.e. 0.1899 0.1773 0.1021 0.1795 0.1735 0.1584
##
## sigma^2 estimated as 186717: log likelihood = -1845.06, aic = 3702.12
#ARIMA(2, 1, 3) gives the minimum AIC value
tsdiag(arima(diff(wholesale_computer_timeseries), order = c(2, 1, 3), include.mean = FALSE))

qqnorm(residuals(arima(diff(wholesale_computer_timeseries), order = c(2, 1, 3), include.mean = FALSE)));qqline(residuals(arima(diff(wholesale_computer_timeseries), order = c(2, 1, 3), include.mean = FALSE)))

ARIMA Forecast
fit_arima = stats::arima(wholesale_computer_timeseries, order = c(2, 1, 3))
plot(forecast(fit_arima, h = 30), pch = 19, ylab = "computer sales")

hold = window(ts(wholesale_computer_timeseries), start = 224)
fit_no_holdout = stats::arima(ts(wholesale_computer_timeseries[-c(224:248)]), order = c(2, 1, 2))
plot(forecast(fit_no_holdout, h = 24))
lines(ts(wholesale_computer_timeseries))

Seasonal ARIMA
plot(window(wholesale_computer_timeseries), start = c(1997,1), ylab = "computer sales")
## Warning in plot.window(xlim, ylim, log, ...): "start" is not a graphical
## parameter
## Warning in title(main = main, xlab = xlab, ylab = ylab, ...): "start" is
## not a graphical parameter
## Warning in axis(1, ...): "start" is not a graphical parameter
## Warning in axis(2, ...): "start" is not a graphical parameter
## Warning in box(...): "start" is not a graphical parameter
Month = c("J", "F", "M", "A", "M", "J", "J", "A", "S", "O", "N", "D")
points(window(wholesale_computer_timeseries, start = c(1997, 1)), pch = Month)

auto.arima(wholesale_computer_timeseries)
## Series: wholesale_computer_timeseries
## ARIMA(3,1,2)(1,0,2)[12] with drift
##
## Coefficients:
## Warning in sqrt(diag(x$var.coef)): NaNs produced
## ar1 ar2 ar3 ma1 ma2 sar1 sma1 sma2
## -0.7351 -0.811 4e-04 0.4938 0.6274 -0.4239 0.1498 -0.419
## s.e. NaN NaN NaN NaN NaN NaN NaN NaN
## drift
## 35.7920
## s.e. 11.6506
##
## sigma^2 estimated as 167516: log likelihood=-1833.77
## AIC=3687.54 AICc=3688.47 BIC=3722.63
#An ARIMA(3,1,2)(1,0,2)[12] is given by auto.arima
plot(diff(diff(wholesale_computer_timeseries),lag=12),xlab='Time',
ylab='First and Seasonal Difference of computer sales')

acf(as.vector(diff(diff(wholesale_computer_timeseries),lag=12)))

#model
m1_wholesale_computer=arima(wholesale_computer_timeseries,order=c(3,1,2),seasonal=list(order=c(1,0,2),
period=12))
m1_wholesale_computer
##
## Call:
## arima(x = wholesale_computer_timeseries, order = c(3, 1, 2), seasonal = list(order = c(1,
## 0, 2), period = 12))
##
## Coefficients:
## ar1 ar2 ar3 ma1 ma2 sar1 sma1 sma2
## -0.6589 -0.7500 0.0446 0.4389 0.6091 -0.5710 0.3275 -0.3995
## s.e. 0.1860 0.1644 0.1006 0.1749 0.1264 0.2037 0.2004 0.0726
##
## sigma^2 estimated as 166222: log likelihood = -1837.29, aic = 3690.57
#Residuals
plot(window(rstandard(m1_wholesale_computer)), start=c(1997,1), ylab='Standardized Residuals', type='o')
## Warning in plot.window(xlim, ylim, log, ...): "start" is not a graphical
## parameter
## Warning in title(main = main, xlab = xlab, ylab = ylab, ...): "start" is
## not a graphical parameter
## Warning in axis(1, ...): "start" is not a graphical parameter
## Warning in axis(2, ...): "start" is not a graphical parameter
## Warning in box(...): "start" is not a graphical parameter
abline(h=0)

acf(as.vector(window(rstandard(m1_wholesale_computer),start=c(1997,1))))

qqnorm(window(rstandard(m1_wholesale_computer),start=c(1997,1)))
qqline(window(rstandard(m1_wholesale_computer),start=c(1997,1)))

shapiro.test(rstandard(m1_wholesale_computer))
##
## Shapiro-Wilk normality test
##
## data: rstandard(m1_wholesale_computer)
## W = 0.98647, p-value = 0.0192
Seasonal ARIMA Forecasting
hold = window(ts(wholesale_computer_timeseries), start = 212)
fit_no_holdout = stats::arima(ts(wholesale_computer_timeseries[-c(212:248)]), order = c(3, 1, 2), seasonal = list(order = c(1, 0, 2), period = 12))
plot(forecast(fit_no_holdout, h = 36))
lines(ts(wholesale_computer_timeseries))

plot(m1_wholesale_computer,n1=c(1997,1),n.ahead=24,xlab='Year',type = 'o',
ylab='Computer Sales')

fit_seasonal_ARIMA = forecast(stats::arima(wholesale_computer_timeseries,order=c(3,1,2),seasonal=list(order=c(1,0,2),period=12)))
fit_seasonal_ARIMA$mean
## Jan Feb Mar Apr May Jun Jul
## 2017
## 2018 20759.58 20853.65 20535.43 20532.17 20321.25 20333.16 20262.69
## 2019 20075.59 20002.97 20093.43 20015.98 20013.27 20064.48 19854.82
## Aug Sep Oct Nov Dec
## 2017 20342.95 20485.66 20868.22 20495.46
## 2018 20249.26 20151.69 20213.96 20142.66 20143.02
## 2019 19583.76
fit_seasonal_ARIMA$lower
## 80% 95%
## Sep 2017 19820.45 19543.86
## Oct 2017 19823.03 19472.26
## Nov 2017 20089.20 19676.81
## Dec 2017 19560.05 19064.87
## Jan 2018 19725.94 19178.76
## Feb 2018 19741.60 19152.91
## Mar 2018 19318.99 18675.04
## Apr 2018 19232.66 18544.74
## May 2018 18956.36 18233.84
## Jun 2018 18888.26 18123.38
## Jul 2018 18744.81 17941.29
## Aug 2018 18672.56 17837.91
## Sep 2018 18539.98 17686.80
## Oct 2018 18560.87 17685.77
## Nov 2018 18455.82 17562.86
## Dec 2018 18424.69 17515.06
## Jan 2019 18320.29 17391.10
## Feb 2019 18214.53 17267.78
## Mar 2019 18274.89 17312.21
## Apr 2019 18163.82 17183.34
## May 2019 18129.04 17131.59
## Jun 2019 18151.11 17138.23
## Jul 2019 17910.36 16881.02
## Aug 2019 17608.53 16562.91
fit_seasonal_ARIMA$upper
## 80% 95%
## Sep 2017 20865.45 21142.04
## Oct 2017 21148.28 21499.06
## Nov 2017 21647.25 22059.64
## Dec 2017 21430.87 21926.05
## Jan 2018 21793.23 22340.40
## Feb 2018 21965.70 22554.39
## Mar 2018 21751.87 22395.81
## Apr 2018 21831.69 22519.61
## May 2018 21686.13 22408.65
## Jun 2018 21778.06 22542.95
## Jul 2018 21780.57 22584.09
## Aug 2018 21825.96 22660.61
## Sep 2018 21763.40 22616.58
## Oct 2018 21867.06 22742.15
## Nov 2018 21829.50 22722.46
## Dec 2018 21861.36 22770.99
## Jan 2019 21830.88 22760.07
## Feb 2019 21791.41 22738.16
## Mar 2019 21911.98 22874.66
## Apr 2019 21868.14 22848.61
## May 2019 21897.51 22894.96
## Jun 2019 21977.84 22990.72
## Jul 2019 21799.28 22828.62
## Aug 2019 21558.98 22604.60
plot(fit_seasonal_ARIMA$residuals)
abline(h = 0)

ARCH/GARCH
diff_wholesales_computer_timeseries = diff(wholesale_computer_timeseries)
plot(diff_wholesales_computer_timeseries);abline(h = 0)

acf(diff_wholesales_computer_timeseries)

pacf(diff_wholesales_computer_timeseries)

acf(abs(diff_wholesales_computer_timeseries))

pacf(abs(diff_wholesales_computer_timeseries))

acf(diff_wholesales_computer_timeseries^2)

pacf(diff_wholesales_computer_timeseries^2)

McLeod.Li.test(y = diff_wholesales_computer_timeseries)

garch_wholesale<-ugarchspec(variance.model = list(model="sGARCH", garchOrder=c(1,1), mean.model = list(armaOrder=c(0,0))), distribution.model = "std")
## Warning: unidentified option(s) in variance.model:
## mean.model
wholesale_garch<-ugarchfit(spec=garch_wholesale, data=diff_wholesales_computer_timeseries)
wholesale_predict<-ugarchboot(wholesale_garch,n.ahead=10, method=c("Partial","Full")[1])
plot(wholesale_predict,which=2)

ga = garch(diff_wholesales_computer_timeseries, order = c(1,1))
##
## ***** ESTIMATION WITH ANALYTICAL GRADIENT *****
##
##
## I INITIAL X(I) D(I)
##
## 1 1.929208e+05 1.000e+00
## 2 5.000000e-02 1.000e+00
## 3 5.000000e-02 1.000e+00
##
## IT NF F RELDF PRELDF RELDX STPPAR D*STEP NPRELDF
## 0 1 1.633e+03
## 1 4 1.633e+03 7.55e-05 2.44e-04 1.1e-07 2.1e+02 4.3e-02 2.57e-02
## 2 5 1.633e+03 2.88e-05 7.54e-05 1.1e-07 2.0e+00 4.3e-02 2.81e-03
## 3 7 1.633e+03 5.97e-06 6.38e-06 8.9e-09 4.5e+00 4.3e-03 6.24e-05
## 4 8 1.633e+03 3.45e-06 4.06e-06 2.2e-08 1.8e+00 8.7e-03 6.86e-06
## 5 10 1.633e+03 1.51e-07 2.77e-07 5.1e-09 5.6e+00 2.2e-03 3.16e-06
## 6 11 1.633e+03 2.55e-08 1.41e-07 5.7e-09 1.3e+00 2.2e-03 1.87e-07
## 7 12 1.633e+03 3.19e-08 3.19e-08 2.1e-09 0.0e+00 8.3e-04 3.19e-08
##
## ***** X-CONVERGENCE *****
##
## FUNCTION 1.632650e+03 RELDX 2.139e-09
## FUNC. EVALS 12 GRAD. EVALS 8
## PRELDF 3.193e-08 NPRELDF 3.193e-08
##
## I FINAL X(I) D(I) G(I)
##
## 1 1.929208e+05 1.000e+00 6.439e-06
## 2 9.873919e-02 1.000e+00 8.839e-04
## 3 1.944549e-02 1.000e+00 -1.439e-04
## Warning in garch(diff_wholesales_computer_timeseries, order = c(1, 1)):
## singular information
qqnorm(ga$residuals);qqline(ga$residuals)
